CD ROM Paradise Collection 4/CD ROM Paradise Collection 4 1995 Nov.iso/science/normix21.zip/NORMIX21/PCNORMIX.FOR

Text File | 1995-01-11 | 75KB | 2,393 lines

c     PC-NORMIX Version 2.1 (January 1995)
c     Copyright (c) 1995 by John H. Wolfe. All rights reserved.
c     This is a revision of Version 1.1 of PCNORMIX for the
c     286 and above. It runs on 386 and 486 pcs, compiled with Microsoft
c     Fortran Power-Station, unix machines with the f77 compiler,
c     and the IBM-4381 mainframe with the FORTVS compiler.
c     Only the I/O differs from the 1978 IBM Mainframe version. The
c     computational algorithms remain unchanged. This version of the program
c     is copyrighted, unlike Version 1.1, which is in the public domain.
c     It may be freely redistributed in its entirety provided that this
c     copyright notice is not removed. It may not be sold for profit or
c     incorporated in commercial programs without the written permission
c     of the copyright holder. Permission is expressly granted for this
c     program to be made available for file transfer from installations
c     offering unrestricted anonymous file transfer on the Internet.
c     This program is provided without any express or implied warranty.
c     The user may modify or recompile the program for his/her own use on any
c     one computer, provided that he/she becomes a Registered user. Distribution
c     of modified or recompiled versions of this program is strictly forbidden
c     without written permission from the copyright holder.
c
c     This (main) program dimensions the arrays.
      implicit real*8 (a-h,o-z)
      real*8 likly
      character*60 fmt
      character*40 nform,dblank
      logical normap
c     external normex
      common/a/likly,xlimit,nvbles,nsampl,ntypes,nsemi,nparam,normap,kx
      common/b/convrg,critrn,ratio,htest,oldlik,en,bigest,thresh,guard,
     1enmin,title(36),fmt,mxhier,maxtyp,maxav,nclust,kscrat,inp,kpun,kv,
     2 lot,lack,luck,iprint,ltime,maxt,iter,iterx,ltye,itmc,nrank,ivalue
     3,method,incase,ntyp(20)
c
c*********************************************************************
c  VAX-UNIX or IBM Mainframe VMS version: Remove C from Column 1
c  on next two lines:
c      dimension a(56000)
c      data nadim/56000/
c     These dimensions are currently set to fit into 571K byte RAM with
c     123K of program and 56K 8-byte words.
c  End of Mainframe statement.
c*********************************************************************
c
c*********************************************************************
c  The next statement works only for MS-Fortran on a PC:
c      dimension a[allocatable,huge](:)
c  The next statement works only for MS-Fortran PowerStation on a PC:
       dimension a[allocatable](:)
c  End of MS-Fortran PC statement.
c*********************************************************************
c
c
c    Logical I/O unit numbers are assigned to files as follows:
c    15      Input-form  statements (filename specified interactively,
c                                    default: "thisjob")
c    11      Input data (filename specified on input-form; if blank,
c                       data are read from unit 15 following the format
c                       statement on the input-form)
c    12      Printout of the analysis (filename specified on input-form,
c                                      default: "prinout")
c     3      "kc3temp" Scratch unit containing factor scores
c     9      "kv9temp" Scratch unit containing raw data
c     4      "discrim" = the discriminant scores
c     7      "dumpout"  the parameter estimates if iteration limit is
c                       reached. These can be used to continue the
c                       analysis by including them as initial estimates
c                       on a new input-form.

      write(*,10)
   10 format(' Enter the file name for the input-form statements (defaul
     1t = "thisjob")')
      read(*,'(a40)',end=20,err=20)nform
      nform = dblank(nform)
      if(nform .eq. ' ')nform='thisjob'
      go to 30
   20 nform = 'thisjob'
   30 open( 15,file=nform)
      write(*,*)' '
      write(*,*)' Reading input form control statements from'
      write(*,*)'          ',nform
      call inform
c
c     calculate number of bytes occupied by each array
      nb=max0(maxav*4,nvbles*8,168)
      nb=((nb+7)/8)*8
      nc=maxtyp*8
c
      nx=nvbles*nsampl*4
      nx=max0(nx,24000,maxav*(maxav+1)*2 )
      nx=((nx+7)/8)*8
c
      nvc=nvbles*nc
      nv=nvbles*8
      np=max0(nsemi*8,nparam*nc)
      if(normap) then
         ncov=nsemi*24
         nvt=nvbles*nv
         nr=nv
         nvv=nvc
      else
         ncov=nsemi*nc
         nvt=8
         nr=8
         nvv=nv
      endif
      kor=2*(nc+nvc+nvt)+ncov+np+nvv+nvbles*nv+(12+nvbles)*nb+ 24
     1 +nsampl*4 + nsampl*24
      kore=kor + nx + 123000
      kor=kor/8 + 1
      korx = nx/8 + 1
      kor = kor + korx
c*********************************************************************
c  The following two statements work only with MS Fortran on a PC.
c  Remove them for the mainframe or unix versions.
      nadim = kor
      allocate (a(kor))
c  End of MS-Fortran PC statements.
c*********************************************************************
      write(*,200) kore,kor,nadim
  200 format(' Storage requirement=',i8,' decimal bytes',/
     1 ' Required:       DIMENSION A( ',i6,')',/
     2 ' This Program:   DIMENSION A( ',i6,')')
      if(kor .gt. nadim)  then
         write(*,*) ' This problem is too large. Execution terminated.'
         stop
      endif
      vbles=nvbles
c     estimate how long this run will take.
      etime=dble(maxav)*(2.865d0*en*vbles+4.961d0*dble(maxav))
      do 40 l=1,maxt
      t=ntyp(l)+1
      iters=min0(iterx,   (10+nvbles/2)*( ntyp(l)-1)+2 )
      ters=iters
      if(normap) then
      etime=etime+16667.d0*dmax1(0.d0,t+t-7.d0)+1.23d0*en*ters*t*
     1(vbles+2.d0)+4.25d0*t*vbles**2+2.12d0*t*vbles**3+67.d0*t*en
      else
      etime=etime+0.62d0*ters*en*t*( vbles+2.0d0)*(vbles+1.d0)+4.25d0*t*
     1vbles**2 +67.d0*en*t
      endif
   40 continue
      etime=etime/25.0d+6
      n1=1
      n2=n1+nc/8
      n3=n2+nvc/8
      n4=n3+ncov/8
      n5=n4+np/8
      n6=n5+nc/8
      n7=n6+nvc/8
      n8=n7+nvv/8
      n9=n8+nvbles**2
      n10=n9+nvt/8
      n11=n10+nvt/8
      n12=n11+nb/8
      n13=n12+nb/8
      n14=n13+nb/8
      n15=n14+nb/8
      n16=n15+(maxav*nvbles+1)/2
      n17=n16+(maxav+1)/2
      n18=n17+nb/8
      n19=n18+nx/8
      n20 = n19 + (nsampl+1)/2
      write (*,300)etime
  300 format(//' Anticipated execution time (on 486/33DX) =',f8.2,
     1' minutes' )
c
      call normex(a,a(n2),a(n3),a(n4),a(n5),a(n6),a(n7),a(n8),a(n9)
     1,a(n10),a(n11),a(n12),a(n13),a(n14),a(n15),a(n16),a(n17),a(n18)
     2,a(n19),a(n20)  )
c
      kc = 3
      close(kc,status='DELETE')
      close(kv,status='DELETE')
      write(*,*)' Analysis completed.'
      stop
c
      end
      subroutine inform
c     This subroutine reads the input form cards on which the user
c     specifies the options controlling the program.
      implicit real*8 (a-h,o-z)
      real*8 likly
      character*60 fmt
      character*40 infile,lsting,dblank
      logical normap
      common/a/likly,xlimit,nvbles,nsampl,ntypes,nsemi,nparam,normap,kx
      common/b/convrg,critrn,ratio,htest,oldlik,en,bigest,thresh,guard,
     1enmin,title(36),fmt,mxhier,maxtyp,maxav,nclust,kscrat,inp,kpun,kv,
     2 lot,lack,luck,iprint,ltime,maxt,iter,iterx,ltye,itmc,nrank,ivalue
     3,method,incase,ntyp(20)
      kpun=7
      maxt=0
      read ( 15, 6) title
6     format ( 9a8)
      read ( 15, 7)nvbles,nsampl,(ntyp(l),l=1,14),nor,enmin,
     1 htest,maxav,mxhier,iterx,iprint,infile,lsting,fmt
  7     format(t21,i2,t37,i4,t65,/t29,14(1x,i2)/t42,i1,t66,f3.0/
     1  t68,f4.3/t34,i3,t65,i3/t20,i3,t42,i1/t22,a40/t20,a40/t13,a60)
      incase = ialong(fmt)
      ltime=600
      if(htest .eq. 0.0d0) htest=0.999d0
      normap=(nor.ne.1)
      if(enmin.eq. 0.0d0) enmin=nvbles+1
      if (iterx.eq.0) iterx=300
      do 2 i=1,14
      if (ntyp(i).ne.0) go to 2
      do 1 j=i,14
      if (ntyp(j).eq.0) go to 1
      ntyp(i)=ntyp(j)
      ntyp(j)=0
      go to 2
1     continue
      maxt=i-1
      go to 3
2     continue
      maxt=14
3     if (maxt.ne.0) go to 5
      maxt=20
      do 4 i=1,20
4     ntyp(i)=i
5     maxtyp=ntyp(maxt)+1
      nsemi=nvbles*(nvbles+1)/2
      nparam=(nvbles+1)*(nvbles+2)/2
      if(normap) nparam=nvbles+1
      en=nsampl
      if(maxav.eq.0) maxav=2000
      maxav=min0(maxav,nsampl)
      ivalue=-1
      infile = dblank(infile)
      if(infile .eq. ' ') then
         inp = 15
      else
         inp = 11
         open( 11,file=infile)
      endif
      lsting = dblank(lsting)
      if(lsting .eq. ' ') lsting = 'prinout'
      call openit(12,lsting)
      write(*,*)' Data will be read on unit#',inp,' from'
      write(*,*)'          ',infile
      write(*,*)' The printout of the analysis will appear on'
      write(*,*)'          ',lsting
      return
      end
c
      function ialong(fmt)
c     Returns ialong = length of the first A format variable specified
c     in the variable format fmt. If ialong = 0, no a-format was found.
c     If ialong = -1, the length of the A format was unspecified.
      character fmt*60,aind*3
      i1 = index(fmt ,'a')
      if(i1 .eq. 0) i1 = index(fmt,'A')
      if(i1 .ne. 0) then
         i2 = index(fmt(i1:i1+3),',')
         leng = i2 - 2
         if (leng .le. 0) then
            print *, 'No length specified for A format'
            ialong = -1
            return
         endif
         aind = fmt(i1+1: i1+leng)
         read(unit=aind,fmt='(i2)')ialong
      else
         ialong = 0
      endif
      return
      end
c

      function dblank(ch)
c     Strip off leading blanks from a character string.
      character*40 ch,dblank,db
      character*1 one(40)
      equivalence (db,one)
      db = ch
      len = 40
      do 10 i=1,len
      if(one(i) .eq. ' ') goto 10
         nb = i - 1
         go to 20
   10 continue
   20 if(nb .eq. 0) go to 50
      do 30 i = 1,len-nb
   30 one(i) = one(i+nb)
      do 40 i = len-nb+1,len
   40 one(i) = ' '
   50 dblank = db
      return
      end
      subroutine normex(lambda,mean,covin,mu,olam,omean,v,t,w,bmat,r,s
     1,b,g,av,dense,member,x,kpath,case)
      implicit real*8 (a-h,o-z)
      real*8 lambda,mean,mu,likly,lnlike(20)
      real*4 x,dense,av
      character*60 fmt
      character*40 fname
      logical normap
      common/a/likly,xlimit,nvbles,nsampl,ntypes,nsemi,nparam,normap,kx
      common/b/convrg,critrn,ratio,htest,oldlik,en,bigest,thresh,guard,
     1enmin,title(36),fmt,mxhier,maxtyp,maxav,nclust,kscrat,inp,kpun,kv,
     2 lot,lack,luck,iprint,ltime,maxt,iter,iterx,ltye,itmc,nrank,ivalue
     3,method,incase,ntyp(20)
      dimension lambda(1), mean(1), covin(1),r(1), s(1), b(1),ndf(20),
     1g(1), mu(1), av(1), olam(1), omean(1), dense(1), x(nvbles,nsampl)
     2,v(1),t(1),w(1),bmat(1),member(1),kpath(2),chi(20),pchi(20)
      logical*1 case(24,nsampl)
c
c     This is the controlling program which calls the other subroutines.
c
c     Notice that, as presently programmed, all data must be in core and
c     stored in the array x. It would be relatively easy to reprogram to
c     handle much larger problems by reading x one record at a time from
c     disk each time initle, histog, moment, and assign are called. One
c     should be careful about this, since x is used as a working storage
c     area in kmean and prplot. Also, this method will slow the program
c     by about 20 percent.
c
      convrg=1.0d-4
      kscrat=4
      fname = 'discrim'
      call openit(kscrat,fname)
      rewind kscrat
      call itime (itma)
c     itime is a routine which returns the clock value in 100-ths of a
c     second.
c     ltime is the time limit in minutes.
      itma=6000*ltime+itma
      write(*,*)' Reading the data'
      call initle (lambda,mean,covin,x,case)
      write(*,*)' Entering HISTOG hierarchical clustering'
      call histog(mean,av,dense,x,r,s,g,b,mu,covin,t,member,kpath,case)
      call itime(itmb)
      tmb=itmb/100.0
c     write (12,15)tmb
      do 8 l=1,maxt
      ltye=l
      iter=0
      method=10 +nvbles
      ntypes=ntyp(l)+1
      write(*,*)' Start REVALU to use new initial estimates'
            call revalu (lambda,mean,covin,av,dense,member,r,s)
      call output (lambda,mean,covin,r,s)
      oldlik=-1.0d08
      call tytle
      ldiv=0
      call itime(itme)
      do 5 iter=1,iterx
c     write(*,*)' Start iteration',iter
            call aitken (lambda,mean,covin,x,r,s,b,v,g,mu,olam,omean)
c      write(*,*)'finish aitken'
      elik=likly-oldlik
      dlik=dabs(elik)
      if (elik.gt.-1.0d-5) go to 1
      ldiv=ldiv+1
      if (ldiv.lt.9) go to 2
c     If the likelihood fails 9 times to exceed the maximum likelihood
c     of prior iterations, terminate job.
      write (12,10)
      if (dlik.gt.convrg) itmc=10
      go to 7
1     oldlik=likly
      ldiv=0
      write (*,11) iter,ntyp(l),likly
c     Test whether the convergence criterion is met.
      if (dlik-convrg) 7,7,2
2     if (iprint.eq.0) go to 3
      call output (lambda,mean,covin,r,s)
      go to 4
3     write (12,11) iter,ntyp(l),likly
4     call itime (itmb)
      itmc=itmb-itma
      tmc=itmc/6000. +ltime
c     Stop and dump current parameter estimates if time limit reached.
c     if (itmc.lt.0) go to 5
c     write (12,13) tmc
c     go to 6
5     continue
      go to 7
6     call dumper (lambda,mean,covin,r,s)
7     call itime(itmb)
c     tme=(itmb-itme)/100.0
c     tmr=tmc*60.0
c     write (12,14)tme,tmr
      call output (lambda,mean,covin,r,s)
      if(ntypes.eq.2) go to 20
      nrank=min0(ntypes-2,nvbles,8)
      write(*,*)' Start ASSIGN'
      call assign(lambda,mean,covin,x,r,s,b,v,g,mu,t,w,bmat,case)
      if(normap) then
                 call mapper(x,r,s)
                 write(*,*)' Finish MAPPER plot of discriminant scores'
      endif
      rewind kv
      read(kv) ((x(i,k),i=1,nvbles),k=1,nsampl)
c     if (itmc.ge.0) go to 9
c     Test whether likelihood is significantly greater for more types.
20    ratio=likly-prelik
      prelik=likly
      ntypes=ntyp(l)
      lnlike(l) = likly
      if (l.le.1) go to 8
      ntyold=ntyp(l-1)
      ndegre=(nparam-1)*(ntypes-ntyold)*2
      chisq=2.0*ratio*(en-0.5-nvbles-ntypes/2.0)/en
      pchisq=chisqp(chisq,ndegre)
      ndf(l) = ndegre
      chi(l) = chisq
      pchi(l) = pchisq
      call tytle
      write (12,12) ntypes,ntyold,ratio,ndegre,chisq,pchisq
      if ((htest .lt. 0.998d0) .and. (pchisq .gt. htest)) go to 9
8     continue
9     if(normap) end file kscrat
      call tytle
      write(12,16)
      write(12,17) ntyp(1),lnlike(1)
      if(maxt .gt. 1) then
      do 90 l=2,maxt
         write(12,17) ntyp(l),lnlike(l),chi(l),ndf(l),pchi(l)
   90 continue
      endif
      return
c
c
10    format (/,5x,'Job terminated--iterations are not converging')
11    format (' Iteration',i4, 3x,'Log likelihood of',i3,' types in this
     1 sample =',g20.8)
12    format (' Logarithm of likelihood ratio of',i3,' to',i3,' types=',
     1 e18.8,/,' Pseudo Chi-square with',i3,' degrees of freedom=',f10.2
     2,/,' "Probability" of null hypothesis=',f11.8)
c 13  format (39h Time limit exceeded. Computation time=f6.2,8h minutes)
c 14  format(//40x,'Iteration time =',f12.2,' seconds.'/40x,'Total time
c    1used =',f11.2,' seconds.')
c 15  format(//40x,'Time in hierarchical grouping =',f8.2,' seconds.')
   16 format(/'         Summary of Likelihood Statistics ',/
     * '  Types    Log Likelihood  Chi-square  DF  "Probability"')
   17 format(i5,g20.8,f12.2,i4,f14.4)
      end
      subroutine openit(no,fname)
      character*40 fname
      logical isther
      inquire(file=fname,exist=isther)
      if(isther) then
      open(no,file=fname,status='OLD')
      else
      open(no,file=fname,status='NEW')
      endif
      return
      end
      subroutine openun(no,fname)
      character*40 fname
      logical isther
      inquire(file=fname,exist=isther)
      if(isther) then
      open(no,file=fname,status='OLD',form='UNFORMATTED')
      else
      open(no,file=fname,status='NEW',form='UNFORMATTED')
      endif
      return
      end
      subroutine initle (lambda,mean,covin,x,case)
      implicit real*8 (a-h,o-z)
      real*8 lambda,mean,likly
      real*4 x
      character*60 fmt*60,acase*24
      logical normap
      common/a/likly,xlimit,nvbles,nsampl,ntypes,nsemi,nparam,normap,kx
      common/b/convrg,critrn,ratio,htest,oldlik,en,bigest,thresh,guard,
     1enmin,title(36),fmt,mxhier,maxtyp,maxav,nclust,kscrat,inp,kpun,kv,
     2 lot,lack,luck,iprint,ltime,maxt,iter,iterx,ltye,itmc,nrank,ivalue
     3,method,incase,ntyp(20)
      dimension lambda(1),mean(nvbles,2),covin(nsemi,2),
     1 x(nvbles,1)
      logical*1 case(24,nsampl),lcase(24)
      equivalence (acase,lcase)
      ntypes=0
      iter=0
      likly=0.0
      do 1 l=1,maxtyp
      lambda(l)=0.0
      do 1 i=1,nvbles
1     mean(i,l)=0.
      naxtyp=maxtyp
      if(normap) naxtyp=3
      do 12 l=1,naxtyp
      do 12 ij=1,nsemi
12    covin(ij,l)=0.0d0
3     do 4 j=1,nsampl
c     Read in all the data.
      if(incase .eq. 0) then
         read (inp,fmt) (x(i,j),i=1,nvbles)
      else
         read (inp,fmt) acase,(x(i,j),i=1,nvbles)
c        write(12,fmt) acase,(x(i,j),i=1,nvbles)
         do 41 i=1,incase
  41     case(i,j) = lcase(i)
      endif
4     continue
      do 5 k=1,nsampl
c     Calculate the overall means and covariances for the whole sample.
      ij=0
      do 5 i=1,nvbles
      mean(i,1)=mean(i,1)+x(i,k)
      do 5 j=i,nvbles
      ij=ij+1
5     covin(ij,1)=covin(ij,1)+x(i,k)*x(j,k)
      ij=0
      do 7 i=1,nvbles
      do 6 j=i,nvbles
      ij=ij+1
6     covin(ij,1)=(covin(ij,1)*en-mean(i,1)*mean(j,1))/(en*en)
7     mean(i,1)=mean(i,1)/en
      return
c
c
8     format (f8.4)
9     format (10f8.4)
      end
      subroutine histog(mean,av,dense,x,r,s,g,b,mu,covin,t,member,kpath,
     * case)
c     The clusters which are formed are used later in the program for
c     initial estimates of the parameters.
      implicit real*8 (a-h,o-z)
      real*8 mean,mu,likly
      real*4 x,dense,r,s,av
      character*60 fmt
      character*40 fname
      logical normap,kmax,grprin
      integer b,g
      common/a/likly,xlimit,nvbles,nsampl,ntypes,nsemi,nparam,normap,kx
      common/b/convrg,critrn,ratio,htest,oldlik,en,bigest,thresh,guard,
     1enmin,title(36),fmt,mxhier,maxtyp,maxav,nclust,kscrat,inp,kpun,kv,
     2 lot,lack,luck,iprint,ltime,maxt,iter,iterx,ltye,itmc,nrank,ivalue
     3,method,incase,ntyp(20)
      logical*1 case(24,nsampl)
      dimension mean(nvbles,maxtyp),x(nvbles,nsampl),av(nvbles,maxtyp),
     1 dense(maxtyp),r(nvbles),g(1),mu(1),covin(nsemi,1),s(1)
     2,member(1),b(1),t(nvbles,nvbles),kpath(1)
      do 2 i=1,nsemi
2     mu(i)=covin(i,1)
      kv=9
      kc=3
      fname = 'kv9temp'
      call openun(kv,fname)
      fname = 'kc3temp'
      call openun(kc,fname)
      do 50 ngrp = 1,3
c     Repeat the hierarchical grouping three times.
      rewind kv
      rewind kc
      write(*,*)'       Factoring the variables'
      call uplow(mu,nvbles,-1)
      call eigen(mu,t,nvbles,0)
      write(kv) ((x(i,k),i=1,nvbles),k=1,nsampl)
      ii=0
      do 4 i=1,nvbles
      ii=ii+i
      ra=1.0d0/dsqrt(mu(ii))
      do 3 j=1,nvbles
3     t(j,i)=t(j,i)*ra
4     continue
      do 21 k=1,nsampl
      do 20 i=1,nvbles
      r(i)=0.0
      do 20 j=1,nvbles
20    r(i)=r(i)+t(j,i)*(x(j,k)-mean(j,1))
c     write(kc) r
      write(kc) (r(ir),ir=1,nvbles)
   21 continue
      write(*,*)'       Ward hierarchical grouping'
      maxhi = mxhier
      if(mxhier .eq. 0) maxhi=30
      maxhi = min0(maxav,maxhi)
      grprin = ((mxhier .ne. 0) .or. (ngrp .eq. 3))
      if(grprin) call tytle
      call kmean(av,dense,x,s,member,b,g,kpath,maxav,maxhi,nsampl,
     1 nvbles,kc,incase,case,grprin)
      write(*,*)'       Compute means of merged clusters'
      rewind kv
      read(kv) ((x(i,k),i=1,nvbles),k=1,nsampl)
      do 24 l=1,maxav
      do 24 i=1,nvbles
24    av(i,l)=0.0
      l=0
      do 28 k=1,nsampl
      ki=kpath(k)
      if(ki .le. 0) then
         ki=-ki
         l=l+1
      endif
      do 28 i=1,nvbles
28    av(i,l)=av(i,l)+x(i,ki)

      ltypes=min0(10,maxav)
      d=0.0
      ij=0
      do 51 i=1,nvbles
      r(i)=0.0
      do 51 j=i,nvbles
      ij=ij+1
   51 mu(ij)=covin(ij,1)+mean(i,1)*mean(j,1)
      kmax=.false.
      do 57 k=1,maxav
      if(member(k) .ge.ltypes) go to 55
      if(d.le. 0.5) go to 55
   58 ij=0
      do 53 i=1,nvbles
      do 53 j=i,nvbles
      ij=ij+1
   53 mu(ij)=mu(ij)-r(i)*r(j)/(d*en)
      if(kmax) go to 50
      d=0.0
      do 54 i=1,nvbles
   54 r(i)=0.0
   55 d=d+dense(k)
      do 56 i=1,nvbles
   56 r(i)=r(i)+av(i,k)
      kmax=k.eq.maxav
      if(kmax) go to 58
   57 continue
   50 continue

      do 40 l=1,maxav
      do 40 i=1,nvbles
   40 av(i,l)=av(i,l)/dense(l)
c      if mxhier = 0, no printout occurs of means at each stage of the
c      hierarchical grouping. Elsewhere in the program, mxhier defaults
c      to 30, but if you want this printout here, you must enter a
c      positive integer explicitly for mxhier on the input form.
      if(mxhier .eq. 0) return
      call tytle
      write (12,1140)
      last=0
      do 440 i=1,maxav
         k=last+1
         l=dense(i)
         last=l+k-1
         write (12,1130)i,l
         write (12,1135)(kpath(j),j=k,last)
440   write (12,1150)(av(j,i),j=1,nvbles)
      do 35 mtypes=2,maxav
         ltypes=maxav-mtypes+2
         last=g(ltypes)
         j=last
         jind=1
         do 29 i=1,nvbles
29       r(i)=av(i,j)*dense(j)
         awt=dense(j)
30       j=j+jind
         if(j.gt.maxav) go to 32
         if(member(j).lt.ltypes*jind) go to 32
         do 31 i=1,nvbles
31       r(i)=r(i)+av(i,j)*dense(j)
         awt=awt+dense(j)
         if(member(j).ge.ltypes) go to 30
32       if(j.lt.last) go to 33
         j=last
         jind=-1
         go to 30
33       do 34 i=1,nvbles
34       r(i)=r(i)/awt
         write (12,1160)last,j
         nj=awt
         nj1=b(j)
         nj2=nj1+nj-1
         write (12,1130)j,nj
         write (12,1135)(kpath(m),m=nj1,nj2)
35    write (12,1150)(r(i),i=1,nvbles)
1130  format(/30x,'Cluster',i4,10x,i4,' points')
1135  format(20i4)
1140  format(42x,'Kmeans analysis'//)
1150  format(10f12.4)
1160  format(/20x,'Merging cluster',i4,' with cluster',i4,
     1' results in the revised' )
      return
c
      end
      subroutine kmean(av,wt,a,best,member,kgrp,lgrp,kpath,n,mxhier,
     1 ns,nv,nti,incase,case,grprin)
c     Adapted from ward's hierarchical grouping program
c     by John H. Wolfe
c
c     This routine combines MacQueen's Kmeans clustering procedure with
c     Ward's hierarchical grouping.
c            ***references***
c     MacQueen,J.   Fifth Berkeley Symposium(1965)
c     Ward,J.  PERSUB Reference Manual(PRL-TR-67-3) (1967)
c     Definitions:
c        n = maximum number of groups
c        ns= number of points to be grouped.
c        nv= number of variables.
c        nti=logical unit containing data to be clustered.
c        wt(i)=number of points in cluster i.
c        av(j,i)=mean of variables j in cluster i.
c        a =  n*(n+1)/2 words of working storage used by the program
c
c         On return from this routine, kpath will be a list of points.
c     kpath(i) will be negative for the first member of one of the n
c     groups and kpath(i+1),kpath(i+2), etc., will have the remaining
c     members in that group.  member(j) will be the number of groups
c     left after the j-th group is merged with another group.
c
c     This routine expects that logical unit nti contains the data to
c     be clustered.  The data should have been properly scaled or
c     converted to factor scores  and written as unformatted (binary)
c     records
c
      logical*1 case(24,ns)
      logical grprin
      dimension a( 1),kgrp(n),kpath(ns),wt(n),member(n),best(n),av(nv,n)
     1,lgrp(n)
c
c     The program first performs a kmeans analysis and then
c     hierarchically groups the groups.
c  1. Initially, each of the first n data points is a separate cluster.
c  2. The inter-group distances, a, are computed.
c  3. The nearest two groups are combined.
c  4. A new data point is read into the space left by the merger.
c  5. The distances of all the clusters to this new point are computed.
c  6. Go back to step 3 until all the data has been read.
c  7. Hierarchically  group the n groups.
c
c     Initialize
      na = (n * (n+1))/2
      ntwo=n+n
      none=n-1
      if(nti.ne.0) rewind nti
      ks=n
      numg=n
      do 10  i=1,n
      if(nti.ne.0) read (nti) (av(k,i),k=1,nv)
c      write(*,*)i,(av(k,i),k=1,nv)
      lgrp(i)=i
      kgrp(i)=i+1
      best(i)=1.0e35
10    wt(i)=1.0
      kgrp(n)=0
      do 15 i=1,ns
 15   kpath(i)=0
c     A group is identified by the least index of its members
c     If i is a member of a group, kpath(i) is another member.
c     If i is the beginning of a group, kgrp(i) begins the next group.
c     lgrp(i)= the leading member of group i
c     A diagonal element of a is the within-group dispersion.
c     The value of a in row i,column j is the dispersion which would
c     result from combining groups i and j .
c     Compute di  block 2
c     write(*,*) 'Compute a(ij)'
      ij=0
      do 40 i=1,none
      ij=ij+1
      a(ij)=0.0
c     write(*,*)i,(av(k,i),k=1,nv)
      do 30 j=i,none
      ij=ij+1
      a(ij)=0.0
      do 20 k=1,nv
 20   a(ij)=a(ij)+(av(k,i)-av(k,j+1))**2
      a(ij)=a(ij)/2.0
      dit=a(ij)
      if(dit.ge.best(i)) go to 30
      best(i)=dit
      member(i)=j+1
 30   continue
 40   continue
      a(na)=0.0
 50   i=1
      di=1.0e35
      do 60 j=1,numg
      if(best(i).ge.di) go to 60
      di=best(i)
      ids=i
 60   i=kgrp(i)
 70   jds=member(ids)
      if(jds.gt.ids) go to 100
      ids=jds
      go to 70
c     Update kgrp and kpath  block 3
100   i=ids
      numg=numg-1
      best(jds)=di
110   if(kgrp(i)-jds)120,130,120
120   i=kgrp(i)
      go to 110
130   kgrp(i)=kgrp(jds)
      kgrp(jds)=-numg
      i=lgrp(ids)
140   if(kpath(i))150,160,150
150   i=kpath(i)
      go to 140
160   kpath(i)=lgrp(jds)
c     Update hyp matrix block 4
      i=1
      best(ids)=1.0e35
      l =((ntwo-ids)*(ids-1))/2+ids
      k =((ntwo-jds)*(jds-1))/2+jds
      m=l+jds-ids
      plop=a(m)
      tsum1=a(m)*(wt(ids)+wt(jds))
      tsum3= (a(l)* wt(ids))+ (a(k)*wt(jds))
      a(l)=plop
      do 350 im=1,numg
      ii=((ntwo-i  )*(i  -1))/2+i
      if(i-ids)210,350,240
210   m=ii+ids-i
      go to 250
240   m=l+i-ids
250   if(i-jds)260,270,270
260   m1=ii+jds-i
      go to 280
270   m1=k+i-jds
280   a(m)=(a(m1)*(wt(i)+wt(jds))+ (
     1a(m)*(wt(i)+wt(ids)))+  (tsum1
     2-tsum3-(a(ii)*wt(i))))/(wt(i)+wt(jds)+wt(ids) )
      dit=a(m)-plop-a(ii)
      if(dit.ge.best(ids)) go to 290
      best(ids)=dit
      member(ids)=i
290   if(member(i).ne.ids.and.member(i).ne.jds) go to 340
      best(i)=1.0e35
      j=i
300   j=kgrp(j)
      if(j.le.0) go to 350
      jj=((ntwo-j)*(j-1))/2+j
      ij=ii+j-i
      dot=a(ij)-a(ii)-a(jj)
      if(dot.ge.best(i)) go to 300
      best(i)=dot
      member(i)=j
      go to 300
340   if(dit.ge.best(i)) go to 350
      best(i)=dit
      member(i)=ids
350   i=kgrp(i)
      wts=wt(ids)
      wt(ids)=wt(ids)+wt(jds)
      if(ks.eq.ns) go to 359
      ks=ks+1
      numg=n
      lgrp(jds)=ks
      kgrp(jds-1)=jds
      kgrp(jds)=jds+1
      kgrp(n)=0
      do 351 i=1,nv
  351 av(i,ids)=(wts*av(i,ids)+wt(jds)*av(i,jds))/wt(ids)
      read(nti) (av(i,jds),i=1,nv)
      wt(jds)=1.0
      best(jds)=1.0e35
      m=jds
      jadd=n
      jj=1
      do 355 j=1,n
      dit=a(jj)
      if(j.eq.jds) go to 354
      do 352 i=1,nv
  352 dit =dit +(av(i,jds)-av(i,j))**2
      a(m)=dit *wt(j  )/(wt(j  )+1.0)
      dit=a(m)-a(jj)
      if(j.gt. jds) go to 353
      if(dit .ge. best(j)) go to 354
      best(j)=dit
      member(j)=jds
      go to 354
353   if(dit .ge. best(jds)) go to 354
      best(jds)=j
      member(jds)=j
      best(jds)=dit
 354  jadd=jadd-1
      if(j.ge.jds) jadd=1
      jj=jj+n-j+1
  355 m=m+jadd
359   if(numg-1) 360,360,50
c     Out put  block 5
360   la=0
      do 361 i=1,n
      l=lgrp(i)
      if(kpath(l).eq.0) la=i
361   kpath(l)=-kpath(l)
      call relist(kpath)
      k=0
      do 362 i=2,ns
      if(kpath(i).ge.0) go to 362
      k=k+1
      do 363 l=1,n
      if(lgrp(l).eq. kpath(i-1)) go to 364
363   continue
364   member(k)=l
      kpath(i-1)=-kpath(i-1)
      kpath(i)=-kpath(i)
      wt(k)=0.0
362   wt(k)=wt(k)+1.0
      wt(n)=wt(n)+1.0
      if(la.eq.0) go to 371
      kpath(ns)=-kpath(ns)
      member(n)=la
      wt(n)=1.0
      do 370 k=1,n
      l=member(k)
366   if(l-k) 367,370,368
367   l=member(l)
      go to 366
368   do 369 j=1,nv
      temp=av(j,k)
      av(j,k)=av(j,l)
369   av(j,l)=temp
370   continue
371   do 372 l=1,n
      k=member(l)
      j=-kgrp(k)
      member(l)=j
372   lgrp(j+1)=l
c     Member(l)= stage of group l
c     lgrp(j)= group whose stage is j-1
      l=0
      do 373 k=1,ns
      if(kpath(k).gt.0) go to 373
      l=l+1
      kgrp(l)=k
373   continue
c     kgrp(l)= the position in kpath which begins group l
      if(grprin) call hgraph(mxhier,n,ns,kpath,member,incase,case)
      return
      end

      subroutine hgraph(mxhier,n,ns,kpath,member,incase,case)
      integer rank,stage
      character*2 bar,gram
      character char*1,alph*24,fmt8x*180,fmt9x*60
      logical*1 case(24,ns),lchar
      dimension kpath(ns),member(n)
      dimension gram(80),bar(2)
      equivalence (char,lchar)
      data bar/' .',' X'/
c     lineln = number of columns displayed on one page
      lineln = 30
      if(incase .gt. 0) then
         m1 = max0(1,(incase+1)/3)
         m2 = incase-m1+2
         write(fmt8x,1010)m1,m2
 1010 format(
     *'(1h,39x,21hHierarchical Grouping,//60x,16hNumber of Groups,/,
     * 3x,24hRank  Item  Kmeans Stage,',i2,'x,2hID,',i2,'x,i3,i8,5i10)')

c      write(12,*)' fmt8x =',fmt8x
      endif

      do 430 lcol=1,mxhier,lineln
         k2=min0(lcol+lineln-1,mxhier)
         linelg = (k2+1-lcol)
         linelg2 = linelg*2
       if(incase .ne. 0) then
c           Set up format fmt9x for one line, depending on part of printout
            write(fmt9x,1020)incase,linelg
 1020       format('(2i6,2i7,4x,a',i2.2,',2x,',i2,'a2)')
c           write(12,*)' fmt9x =',fmt9x
         endif

         kmean=0
         do 420 rank=1,ns
            if(mod(rank,50) .eq. 1) then
               if(incase .eq. 0) then
                 write (12,1030)lcol,(j,j=lcol+4,k2,5)
              else
                 write(12,fmt8x)lcol,(j,j=lcol+4,k2,5)
               endif
            endif
 1030 format(1h,39x,'Hierarchical Grouping'//60x,'Number of Groups'/3x
     1,'Rank  Item  Kmeans Stage',5x,i10,i8,5i10)

            stage=n+1
            j=kpath(rank)
            if(j.le.0) then
               kmean=kmean+1
               stage=member(kmean)
               j=-j
            endif
            lk=max0(0,stage -lcol+1)+ 1
            if(incase .gt. 0) then
c              Copy case ==> alph
               do 400 ii=1,incase
                  lchar = case(ii,j)
  400          alph(ii:ii) = char
            endif

            do 410 k=1,linelg2
               if(k .lt. lk) then
                  gram(k)=bar(2)
               else
                  gram(k)=bar(1)
               endif
  410       continue

            if(incase .eq. 0) then
               write (12,1040)rank,kpath(rank),kmean,stage,
     *         (gram(ig),ig=1,linelg)
 1040          format(2i6,2i7,14x,30a2)
            else
               write (12,fmt9x)rank,kpath(rank),kmean,stage,alph,
     *         (gram(ig),ig=1,linelg)
            endif
  420    continue
  430 continue


      return
      end

      subroutine relist(kpath)
c
c     Given a list, kpath, such that kpath(i) is the next element of
c     the list after i, this routine reorders kpath so that kpath(i) is
c     the i-th element of the list.
c
      dimension kpath(1)
      l=1
      j=0
1     j=j+1
      m=iabs(l)
2     if(m-j) 3,4,5
3     m=iabs(kpath(m))
      go to 2
4     i=kpath(m)
      go to 6
5     i=kpath(m)
      kpath(m)=kpath(j)
6     kpath(j)=l
      l=i
      if(i.ne.0) go to 1
      return
      end

      subroutine revalu (lambda,mean,covin,av,dense,member,r,s)
c     This subroutine re-evaluates previous estimates of the parameters
c     and finds initial estimates preparatory to iteration. It is called
c     every time the number of hypothesized types is increased.
c     Revised 9 Jan 95 to fix a bug reported by B. Robinson in which
c     only the means of the first cluster were read under the normap option.
      implicit real*8 (a-h,o-z)
      real*8 lambda,mean,likly
      real*4   dense ,av
      character*60 fmt
      logical normap
      common/a/likly,xlimit,nvbles,nsampl,ntypes,nsemi,nparam,normap,kx
      common/b/convrg,critrn,ratio,htest,oldlik,en,bigest,thresh,guard,
     1enmin,title(36),fmt,mxhier,maxtyp,maxav,nclust,kscrat,inp,kpun,kv,
     2 lot,lack,luck,iprint,ltime,maxt,iter,iterx,ltye,itmc,nrank,ivalue
     3,method,incase,ntyp(20)
      dimension lambda(1), mean(nvbles,2), covin(nsemi,2), av(nvbles,2),
     1 dense(1),member(1),r(nvbles),s(nvbles)
      ltypes=ntypes-1
      do 1 l=2,ntypes
      lambda(l)=0.0d0
      do 1 i=1,nvbles
1     mean(i,l)=0.0d0
      l=1
      do 2 k=1,maxav
      if(member(k).lt.ltypes) l=l+1
      lambda(l)=lambda(l)+dense(k)
      do 2 i=1,nvbles
2     mean(i,l)=mean(i,l)+av(i,k)*dense(k)
      do 4 l=2,ntypes
      do 3 i=1,nvbles
3     mean(i,l)=mean(i,l)/lambda(l)
4     lambda(l)=lambda(l)/en
      ij=0
      do 5 i=1,nvbles
      do 5 j=i,nvbles
      ij=ij+1
      covin(ij,2)=covin(ij,1)+mean(i,1)*mean(j,1)
      do 5 l=2,ntypes
5     covin(ij,2)=covin(ij,2)-lambda(l)*mean(i,l)*mean(j,l)
      if(.not. normap) then
         do 6 l=2,ntypes
         do 6 i=1,nsemi
6        covin(i,l)=covin(i,2)
      endif
    8 if((ivalue .eq. 0) .or. (ltypes .lt. ivalue)) return
      if(ltypes .eq. ivalue) go to 10
      read ( 15, 100)ivalue,imn,istd,icor
      go to 8
   10 sump = 0.0d0
      do 17 l=2,ntypes
      read ( 15, *)p
      if(p.ne.0.0d0) lambda(l)=p
      sump = sump + lambda(l)
      if(imn.ne.0) read ( 15, *)(mean(i,l),i=1,nvbles)
      if(normap .and. l .ne. ntypes) go to 17
      if(istd.ne.0) then
         read ( 15, *)(s(i),i=1,nvbles)
         k=0
         do 15 i=1,nvbles
         if(icor.eq.0) then
            do 12 j=1,nvbles
   12       r(j)=0.0d0
            r(i)=1.0d0
         else
            read ( 15, *)(r(j),j=1,nvbles)
       endif
         do 15 j=i,nvbles
         k=k+1
   15    covin(k,l)=r(j)*s(j)*s(i)
      endif
   17 continue
c     Standardize lambdas to sum to one.
      do 18 l=2,ntypes
   18 lambda(l) = lambda(l)/sump
      return
  100 format(25x,i2,16x,i1,10x,i1,14x,i1)
      end

      subroutine itime(it)
      it = 0
      return
      end

      subroutine aitken (lambda,mean,covin,x,r,s,b,v,g,mu,olam,omean)
c     This subroutine does all the computations of a single iteration
c     for improved estimates of the parameters.
      implicit real*8 (a-h,o-z)
      real*8 lambda,mean,mu,likly
      real*4 x
      character*60 fmt
      logical normap
      common/a/likly,xlimit,nvbles,nsampl,ntypes,nsemi,nparam,normap,kx
      common/b/convrg,critrn,ratio,htest,oldlik,en,bigest,thresh,guard,
     1enmin,title(36),fmt,mxhier,maxtyp,maxav,nclust,kscrat,inp,kpun,kv,
     2 lot,lack,luck,iprint,ltime,maxt,iter,iterx,ltye,itmc,nrank,ivalue
     3,method,incase,ntyp(20)
      dimension lambda(1),mean(nvbles,2),covin(nsemi,2),mu(nparam,2),
     1   b(1),olam(1),omean(nvbles,2) ,r(1),s(1),g(1),x(1),v(nvbles,1)
      xlimit=amin0(40,(iter/5)*10)
c     When xlimit=0, all points are assigned a membership probability of
c     1.0 in their nearest cluster.
      n=nvbles-1
      luck=1
      bigest=20.0
      small=1.0-1.0/bigest
      lqq=nsemi
      do 3 l=2,ntypes
      if(l.gt.2.and.normap) go to 2
      call syminv (covin,nvbles,detm,r,s,g,lqq)
      if(normap) go to 2
      if (enmin.lt.en*lambda(l)) go to 2
      lsing=l-1
      nsing=en*lambda(l)
      oldlik=-1.0d+10
      write (12,101)lsing,nsing,oldlik
  101 format(' Type',i3,' with ',i3,' points approaching singularity'/
     1' Log Likelihood re-initialized to',d15.8)
      ij=0
      do 1 i=1,nvbles
      do 1 j=i,nvbles
      ij=ij+1
      covin(ij,l)=covin(ij,1)+mean(i,1)*mean(j,1)
      do 1 m=2,ntypes
1     covin(ij,l)=covin(ij,l)-mu(i+1,m)*mu(j+1,m)/(mu(1,m)*en)
      call syminv (covin,nvbles,detm,r,s,g,lqq)
2       b(l)=dlog(lambda(l))-0.5*detm
3     lqq=lqq+nsemi
      likly=0.0
      do 4 l=2,ntypes
      do 4 i=1,nparam
4     mu(i,l)=0.0
      if(normap) go to 40
      iii=nvbles
      ii=1
      do 6 i=1,nvbles
      do 5 l=2,ntypes
5     covin(ii,l)=covin(ii,l)/2.0
      ii=ii+iii
6     iii=iii-1
      go to 47
40    k=1
      do 42 i=1,nvbles
      likly=likly-0.5*covin(k,2)*(covin(k,1)+mean(i,1)**2)*en
      k=k+1
      if (i.gt.n) go to 42
      do 41 j=i,n
      likly=likly-covin(k,2)*(covin(k,1)+mean(j+1,1)*mean(i,1))*en
41    k=k+1
42    continue
      do 46 l=2,ntypes
      do 45 i=1,nvbles
      k=i
      n=nvbles
      v(i,l)=0.0d0
      if (i.eq.1) go to 44
      do 43 j=2,i
      v(i,l)=v(i,l)+covin(k,2)*mean(j-1,l)
      n=n-1
      k=n+k
43    continue
44    do 45 j=i,nvbles
      v(i,l)=v(i,l)+covin(k,2)*mean(j  ,l)
45    k=k+1
      do 46 i=1,nvbles
46    b(l)=b(l)-v(i,l)*mean(i,l)*0.5d0
c     write(*,*)'Start moment'
47          call moment (mean,covin,x,r,s,b,g,mu,v)
c     write(*,*)'Finish moment'
      byg=1.0
      if (method.gt.0) go to 8
      byg=bigest
      do 7 l=2,ntypes
      delta=mu(1,l)/en-lambda(l)
      biga=lambda(l)-olam(l)
      sml=1.0-1.0/byg
      if (dabs(delta).lt.sml*dabs(biga)) byg=biga/(biga-delta)
      if(delta) 21,7,22
21    bound=dmax1(lambda(l)*0.5d0,1.0d0/en)
      if(bound.lt.-delta*byg) byg=-bound/delta
      go to 7
22    bound=dmin1((1.0d0-lambda(l))*0.5d0,(en-1.0d0)/en)
      if(bound.lt.delta*byg) byg=bound/delta
7     continue
8     do 12 l=2,ntypes
      if(mu(1,l).lt. 1.0d0) mu(1,l)=1.0d0
      olam(l)=lambda(l)
      lambda(l)=lambda(l)+byg*(mu(1,l)/en-lambda(l))
      lambda(l)=dmax1(1.0d-1/en,lambda(l))
      ij=0
      do 9 i=1,nvbles
      biga=mean(i,l)-omean(i,l)
      bigc=mu(i+1,l)/mu(1,l)-mean(i,l)
      big=bigest
      if (dabs(bigc).lt.small*dabs(biga)) big=biga/(biga-bigc)
      if (method.gt.0) big=1.0
      omean(i,l)=mean(i,l)
9     mean(i,l)=mean(i,l)+big*bigc
      if(normap) go to 12
      kj=nvbles+1
      do 10 i=1,nvbles
      do 10 j=i,nvbles
      ij=ij+1
      kj=kj+1
10    covin(ij,l)=(mu(kj,l)-mu(i+1,l)*mean(j,l)-mu(j+1,l)*mean(i,l))/mu(
     11,l)+mean(i,l)*mean(j,l)
12    continue
      method=iabs(method-1)
      if (likly.lt.oldlik) method=luck
      if(.not.normap) return
      do 11 i=1,nvbles
      do 11 j=i,nvbles
      ij=ij+1
      covin(ij,2)=covin(ij,1)+mean(i,1)*mean(j,1)
      do 11 l=2,ntypes
11    covin(ij,2)=covin(ij,2)-mu(i+1,l)*mu(j+1,l)/(mu(1,l)*en)
      return
      end
      subroutine moment (mean,covin,x,r,s,b,g,mu,v)
      implicit real*8 (a-h,o-z)
      real*8 mean,mu,likly
      real*4 x
      logical normap
      common/a/likly,xlimit,nvbles,nsampl,ntypes,nsemi,nparam,normap,kx
      dimension mean(1),covin(1),x(1),r(1),s(1),b(1),g(1),mu(1),v(1)
c     This routine computes
c     mu(1,l)= sum over all individuals of their probabilities of
c     membership in type l.
c     mu(i+1,l)= sum over all individuals of their x(i,k)
c     times their probability of membership in type l.
c     mu(ij,l) = sum over all individuals of their x(i,k)*x(j,k)
c     times their probability of membership in type l.
c     g(l)=probability of membership in type l,
c
c
c     This subroutine contains the most deeply nested loops of the
c     entire normix program.  The following tricks are used
c     to improve object code efficiency:
c     --- The routine computes its own indexes for the
c     doubly-indexed arrays instead of depending on the compiler.
c     x(ik)  is  x(i,k)
c     mean(il)   is mean (i,l)
c     covin(m) is covin(i,j,l)
c     mu(ij)  is mu(i,j,l)
      do 12 k=1,nsampl
            call densit(mean,covin,x,r,s,b,v,g,k)
      if(normap) go to 44
      ij=nparam
      do 11 l=2,ntypes
      pm=g(l)
      ia=ij+1
      ij=ia+nvbles
      mu(ia)=pm+mu(ia)
      do 11 i=1,nvbles
      z=r(i)*pm
      ia=ia+1
      mu(ia)=z+mu(ia)
      do 11 j=i,nvbles
      ij=ij+1
11    mu(ij)=mu(ij)+z*r(j)
      go to 12
44    ia=nparam
      do 46 l=2,ntypes
      pm=g(l)
      ia=ia+1
      mu(ia)=pm+mu(ia)
      do 46 i=1,nvbles
      ia=ia+1
46    mu(ia)=r(i)*pm+mu(ia)
12    continue
      return
      end
      subroutine tytle
      implicit real*8 (a-h,o-z)
      real*8 likly
      character*60 fmt
      logical normap
      common/a/likly,xlimit,nvbles,nsampl,ntypes,nsemi,nparam,normap,kx
      common/b/convrg,critrn,ratio,htest,oldlik,en,bigest,thresh,guard,
     1enmin,title(36),fmt,mxhier,maxtyp,maxav,nclust,kscrat,inp,kpun,kv,
     2 lot,lack,luck,iprint,ltime,maxt,iter,iterx,ltye,itmc,nrank,ivalue
     3,method,incase,ntyp(20)
      write (12,1)
      if(normap) write (12,2)
      write (12,3)title
      return
c
c
1     format (1h,42x,'PC-NORMIX  Version 2.1',/,25x,
     *'Wolfe Normal Mixture Analysis Procedure (January 1995 Revision)')
2     format(31x,'Equal Covariance Matrices Within Groups')
3     format(//4(1x, 9a8/))
      end
      subroutine output (lambda,mean,covin,r,s)
      implicit real*8 (a-h,o-z)
      real*8 lambda,mean,likly
      character*60 fmt
      logical normap
      common/a/likly,xlimit,nvbles,nsampl,ntypes,nsemi,nparam,normap,kx
      common/b/convrg,critrn,ratio,htest,oldlik,en,bigest,thresh,guard,
     1enmin,title(36),fmt,mxhier,maxtyp,maxav,nclust,kscrat,inp,kpun,kv,
     2 lot,lack,luck,iprint,ltime,maxt,iter,iterx,ltye,itmc,nrank,ivalue
     3,method,incase,ntyp(20)
      dimension mean(nvbles,2), covin(nsemi,2), lambda(1), r(1), s(1)
      call tytle
      ntp=ntypes-1
      ktypes=ntypes
      if(normap) ktypes=2
c     write (*,9) nsampl,nvbles,ntp,iter,ntp,likly
      write (12,9) nsampl,nvbles,ntp,iter,ntp,likly
      do 8 l=1,ktypes
      if (l-1) 1,1,2
1     write (12,10)
      go to 3
2     ltyp=l-1
      if(normap) go to 20
      write (12,11) ltyp,lambda(l)
      go to 3
20    write (12,21)
 3    lqq=1
      do 4 i=1,nvbles
      s(i)=dsqrt(covin(lqq,l))
4     lqq=lqq+nvbles-i+1
      if(normap.and.l.ne.1) go to 28
      write (12,12) (i,i=1,nvbles)
      write (12,13)(mean(i,l),i=1,nvbles)
28    write (12,14) (i,i=1,nvbles)
      write (12,13) (s(i),i=1,nvbles)
      write (12,15) (i,i=1,nvbles)
      do 8 i=1,nvbles
      kj=i-1
      ij=i
      if (kj.eq.0) go to 6
      do 5 j=1,kj
      r(j)=covin(ij,l)/(s(i)*s(j))
5     ij=ij+nvbles-j
6     do 7 j=i,nvbles
      r(j)=covin(ij,l)/(s(i)*s(j))
7     ij=ij+1
8     write (12,13) (r(j),j=1,nvbles)
      if(.not. normap) return
      write (12,16) (i,i=1,nvbles)
      do 29 l=2,ntypes
      ltyp=l-1
29    write (12,17) ltyp,lambda(l),(mean(i,l),i=1,nvbles)
      return
c
c
9     format (/10x,'Sample Size =',i11,/,10x,'Number of Variables =',i3,
     1/,10x,'Number of Types =',i7,//25x,'Iteration Number',i4,//3x,
     2'Log Likelihood of',i3,' Types in this sample =',g20.8,/)
10    format (20x,'Characteristics of the Whole Sample',//)
11    format (/20x,'Characteristics of Type',i4,//15x,'The Proportion of
     1 the Population from this Type=',f6.3/)
12    format (34x,'Means',/10(i8,9i12,/))
13    format (10f12.4)
14    format (/27x,'Standard Deviations',/10(i8,9i12,/))
15    format (/32x,'Correlations',/10(i8,9i12,/))
21    format (//20x,'Within-group Standard Deviations and Correlations')
16    format (//30x,'Means',/' Type Lambda',i6,9i10,/10(12x,i6,9i10,/))
17    format (i4,f8.4,10f10.4/10(12x,10f10.4/))
      end
      subroutine dumper (lambda,mean,covin,r,s)
      implicit real*8 (a-h,o-z)
      real*8 lambda,mean,likly
      character*60 fmt
      character*40 fname
      logical normap
      common/a/likly,xlimit,nvbles,nsampl,ntypes,nsemi,nparam,normap,kx
      common/b/convrg,critrn,ratio,htest,oldlik,en,bigest,thresh,guard,
     1enmin,title(36),fmt,mxhier,maxtyp,maxav,nclust,kscrat,inp,kpun,kv,
     2 lot,lack,luck,iprint,ltime,maxt,iter,iterx,ltye,itmc,nrank,ivalue
     3,method,incase,ntyp(20)
      dimension mean(nvbles,2), covin(nsemi,2), lambda(1), r(1), s(1)
      n=2
c     This routine writes parameters to a file in format readable by the
c     subroutine initle.
      fname ='dumpout'
      call openit(kpun,fname)
      do 5 l=2,ntypes
      if(.not. normap) n=l
      ltyp=l-1
      write(kpun,6) lambda(l),ltyp
      write(kpun,7) (mean(i,l),i=1,nvbles)
      lqq=1
      do 1 i=1,nvbles
      s(i)=dsqrt(covin(lqq,n))
1     lqq=lqq+nvbles-i+1
      write(kpun,7) (s(i),i=1,nvbles)
      do 5 i=1,nvbles
      kj=i-1
      ij=i
      if (kj.eq.0) go to 3
      do 2 j=1,kj
      r(j)=covin(ij,n)/(s(i)*s(j))
2     ij=ij+nvbles-j
3     do 4 j=i,nvbles
      r(j)=covin(ij,n)/(s(i)*s(j))
4     ij=ij+1
5     write(kpun,7) (r(j),j=1,nvbles)
      return
c
6     format (f8.4,5x,'type',i3)
7     format (10f8.4)
      end
      subroutine densit (mean,covin,x,r,s,b,v,g,k)
      implicit real*8 (a-h,o-z)
      real*8 mean,likly
      real*4 x
      character*60 fmt
      logical normap
      common/a/likly,xlimit,nvbles,nsampl,ntypes,nsemi,nparam,normap,kx
      common/b/convrg,critrn,ratio,htest,oldlik,en,bigest,thresh,guard,
     1enmin,title(36),fmt,mxhier,maxtyp,maxav,nclust,kscrat,inp,kpun,kv,
     2 lot,lack,luck,iprint,ltime,maxt,iter,iterx,ltye,itmc,nrank,ivalue
     3,method,incase,ntyp(20)
      dimension mean(1),covin(1),x(1),r(1),s(1),b(1),g(2),v(1)
      if(k .gt. 1) go to 50
      xlimit=40.0d0
      n=nvbles-1
      assign 3 to lone
      if (n.eq.0) assign 6 to lone
      ik=1
50    eph=0.0d0
      do 1 i=1,nvbles
      r(i)=x(ik)
1     ik=ik+1
      il=nvbles
      if(normap) go to 40
      ij=nsemi
      do 7 l=2,ntypes
      do 2 i=1,nvbles
      il=il+1
2     s(i)=r(i)-mean(il)
      dsqrd=  b(l)
      go to lone, (3,6)
3     do 5 i=1,n
      ij=ij+1
      cute=s(i)*covin(ij)
      do 4 j=i,n
      ij=ij+1
4     cute=cute+covin(ij)*s(j+1)
5     dsqrd=dsqrd-cute*s(i)
6     ij=ij+1
7     g(l)=dsqrd-covin(ij)*s(nvbles)**2
      go to 43
40    do 42 l=2,ntypes
      dsqrd= b(l)
      do 41 i=1,nvbles
      il=il+1
41    dsqrd=dsqrd+v(il)*r(i)
42    g(l)=dsqrd
43    biga=g(2)
      lot=2
      bigb=biga
      if (ntypes.eq.2) go to 13
      do 9 l=3,ntypes
      if (g(l).le.bigb) go to 9
      if (g(l).le.biga) go to 8
      bigb=biga
      biga=g(l)
      lot=l
      go to 9
8     bigb=g(l)
9     continue
      if ((biga-bigb).gt.xlimit) go to 13
      g(lot)=1.0d0
      do 10 l=2,ntypes
      if (l.eq.lot) go to 10
      dsqrd=g(l)-biga
      g(l)=0.0d0
      if (dsqrd.gt.-4d1   ) g(l)=dexp(dsqrd)
10    eph=eph+g(l)
      likly=likly+dlog(eph)+biga
      do 11 l=2,ntypes
11    g(l)=g(l)/eph
      return
c     Consider the probability of membership to be one for nearest type
13    likly=likly+biga
      do 14 l=2,ntypes
14    g(l)=0.0d0
      g(lot)=1.0d0
      return
      end
c ***********************************************************************
      subroutine assign(lambda,mean,covin,x,r,s,b,v,g,mu,t,w,bmat,case)
c     This routine prints the probabilities of membership of each
c     individual in each type.
c     It computes and prints discriminant functions and each individuals
c     discriminant scores.
c      Discriminant analysis.  See the following references
c     Cooley+Lohnes. Multivariate Proced.Behav.Sciences  pages 117,189.
c     Friedman and Rubin. Journ. Am. Stat. Assn. 1967, page 1167.
c
      implicit real*8 (a-h,o-z)
      real*8 lambda,mean,mu,likly
      real*4 x
      integer is(12)
      character*72 fmtx1,fmtx2,fmtx3,fmtx4,fmtx5,fmtx6,fmtx7
      character char*1,idstrng*3,fmt*60
      character*24 alph,id(12)
      logical*1 case(24,nsampl),lchar
      logical normap,newhead

      common/a/likly,xlimit,nvbles,nsampl,ntypes,nsemi,nparam,normap,kx
      common/b/convrg,critrn,ratio,htest,oldlik,en,bigest,thresh,guard,
     1enmin,title(36),fmt,mxhier,maxtyp,maxav,nclust,kscrat,inp,kpun,kv,
     2 lot,lack,luck,iprint,ltime,maxt,iter,iterx,ltye,itmc,nrank,ivalue
     3,method,incase,ntyp(20)
      dimension covin(nsemi,1),lambda(1),g(1),r(nrank),s(1),mean(1),b(1)
     1,v(1),t(nvbles,nvbles),w(nvbles,nvbles),bmat(nvbles,nvbles),
     1 x(nvbles,1),mu(nparam,1)
      equivalence (char,lchar)
c
c
      call itime(itmda)
      kc=3
      rewind kc
      al=likly
      if(normap) then
        mrank= min0(ntypes-2,nvbles)
         ij=0
         do 10 i=1,nvbles
         do 10 j=i,nvbles
            ij=ij+1
            w(i,j)=covin(ij,2)
            bmat(i,j)=covin(ij,1)-covin(ij,2)
   10   bmat(j,i)=bmat(i,j)
c        bmat is the between-group covariance matrix
c        w is the within-group covariance matrix
         call nroot(nvbles,bmat,w,r,t)
         call tytle
         write (12,3000)(l,l=1,mrank)
 3000   format (//25x,'Discriminant Functions',/(16x,10i10))

         do 20 i=1,nvbles
   20   write (12,3010) i,(t(i,l),l=1,mrank)
 3010    format (' Variable',i4,7x,10f10.4/(20x,10f10.4))

         write (12,3020) (r(l),l=1,mrank)
 3020   format (/' Eigenvalues of B/W=',10f10.4)

         call itime(itmdb)
         tmd=(itmdb-itmda)/100.0
c        write (12,3030)tmd
c3030    format(//40x,'Time to Compute Discriminant Functions =',f8.2,
c    1' seconds.')
      else
         lqq=nsemi
         do 40 l=2,ntypes
         ij=0
         do 30 i=1,nvbles
         do 30 j=i,nvbles
         ij=ij+1
   30    mu(ij,l)=covin(ij,l)
         call syminv (covin,nvbles,detm,r,s,g,lqq)
           b(l)=dlog(lambda(l))-0.5*detm
   40    lqq=lqq+nsemi
         iii=nvbles
         ii=1
         do 60 i=1,nvbles
         do 50 l=2,ntypes
   50   covin(ii,l)=covin(ii,l)/2.0
         ii=ii+iii
   60 iii=iii-1
      endif

      likly=0.0d0
      ntp=ntypes-1
      if(incase .gt. 0) idstrng = ' ID'
      if(incase .eq. 0) idstrng = '   '
      m1 = (100-incase)/6
      lineper = (ntp + m1 -1)/m1
      m3 = 20 + incase
      m4 = 5
      m5 = max0(1,incase/3)
      m5r = max0(1,incase - m5)
      m6 = 120/(13+incase)
      m7 = m5 + m5r

c     Setup format fmtx1 for probability title.
         write(unit=fmtx1,fmt=3040)m5+1,(incase+12)
 3040 format('(8h   Seq.#,',i2,'x,a3,t',i2.2,',7hCluster,8x,',
     * '24hMembership Probabilities)')

c     Setup format fmtx2 for discriminant score title.
         if(normap) write(unit=fmtx2,fmt=3050)m5+1,(incase+12)
 3050 format('(8h   Seq.#,',i2,'x,a3,t',i2.2,',7hCluster,8x,',
     * '19hDiscriminant Scores)')

c     Setup format fmtx3 for cluster numbers at top of table.
         write(unit=fmtx3,fmt=3060)(incase+19),m1
 3060    format( '(t',i2.2,',',i2,'i6)' )
c     Setup format fmtx4 for probabilities or discriminant scores.
         if(incase .eq. 0) then
            fmtx4 = '(i6,i9,4x,16f6.3,:,(/19x,16f6.3,:))'
            fmtx5 = '(i6,i9,4x,16f6.2,:,(/19x,16f6.2,:))'
         else
            write(unit=fmtx4,fmt=3070)incase,m1,m3,m1
 3070       format('(i6,3x,a',i2.2,',i7,4x,',i2,'f6.3,:, (/',i2.2,'x,',
     *               i2,'f6.3,:))')
            write(unit=fmtx5,fmt=3080)incase,m1,m3,m1
 3080       format('(i6,3x,a',i2.2,',i7,4x,',i2,'f6.2,:, (/',i2.2,'x,',
     *               i2,'f6.2,:))')
      endif
c     Set up format for printing case lists sorted by cluster.
      write(fmtx6,4020)m6,m5 ,m5r
 4020 format( '(',i2,'(8h    Seq#,',i2,'x,a3,',i2,'x))' )
      if(incase .eq. 0) then
         write(fmtx7,4030)m6,m7
 4030    format('(',i2,'(i8,3x,',i2,'x))')
      else
         write(fmtx7,4040)m6,incase
 4040    format('(',i2,'(i8,3x,a',i2.2,'))')
      endif


c      write(12,*)' fmtx1=',fmtx1
c      write(12,*)' fmtx2=',fmtx2
c      write(12,*)' fmtx3=',fmtx3
c      write(12,*)' fmtx4=',fmtx4
c      write(12,*)' fmtx5=',fmtx5
c      write(12,*)' fmtx6=',fmtx6
c      write(12,*)' fmtx7=',fmtx7
c *****************************************************************
      do 100 k=1,nsampl
c     if line count > 44, then start a new page.
      if(mod(k,45/lineper) .eq. 1) then
         call tytle
         write (12,fmtx1)idstrng
         write (12,3090)
3090      format( 1h  )
c
         write (12,fmtx3) (i,i=1,ntypes-1)
         write (12,3090)
      endif

      if(incase .gt. 0) then
c     Copy case ==> alph
         do 70 i=1,incase
            lchar = case(i,k)
   70    alph(i:i) = char
      endif

c     Compute the membership probabilities
      call densit (mean,covin,x,r,s,b,v,g,k)
c     Find the largest probability, and assign case.
      xmax=0
      do 80 l=2,ntypes
         if (g(l).lt.xmax) go to 80
         lot=l
c         lot = the most probable type for individual k.
         xmax=g(l)
   80 continue
      lit = lot -1
c     If normap, then compute discriminant scores.
      if(normap) then
         do 90 l=1,nrank
         r(l)=0.0
         do 90 i=1,nvbles
   90    r(l)=r(l)+t(i,l)*(x(i,k)-mean(i))
         write(kc)k,lot,(r(i),i=1,nrank)

c     Save discriminant scores on tape, in card-image format.
         write(kscrat,4000)k,lit,(r(i),i=1,nrank)
 4000    format (i6,i4,10f7.2)
      else
         write(kc)k,lot
      endif

c     Printout probabilities of membership
         if(incase .eq. 0) then
               write(12,fmtx4)k,lit,(g(i),i=2,ntypes)
         else
               write(12,fmtx4)k,alph,lit,(g(i),i=2,ntypes)
         endif

  100 continue
      likly=al
      if(normap) then
c *****************************************************************
c     Print the discriminant scores
      rewind kc
      do 120 k=1,nsampl
c     Start a new page if the line count > 44.
      if(mod(k,45) .eq. 1) then
         call tytle
         write (12,fmtx2)idstrng
         write (12,3090)
c
         write (12,fmtx3) (i,i=1,nrank)
         write (12,3090)
      endif
      if(incase .gt. 0) then
c     Copy case ==> alph
         do 110 i=1,incase
            lchar = case(i,k)
  110    alph(i:i) = char
      endif
c     Read the cluster and discriminant scores from unit kc.
      read(kc)kk,lot,(r(i),i=1,nrank)
      lit = lot -1
c     Print the current case.
      if(incase .eq. 0) then
         write(12,fmtx5)k,lit,(r(i),i=1,nrank)
      else
         write(12,fmtx5)k,alph,lit,(r(i),i=1,nrank)
      endif
  120 continue
c     If not normap, then
      else
         do 130  l=2,ntypes
         do 130 i=1,nsemi
  130    covin(i,l)=mu(i,l)
      endif
c**********************************************************************
c     Print case numbers and IDs for each cluster.
      lap=45
      do 160 lotpas=1,ntp
      ik = 0
      newhead = .true.
      rewind kc
c     Search through all the cases, looking for ones in this cluster.
      do 150 k=1,nsampl
c        Print new page if line count > 44.
      if(lap .ge. 45) then
            lap = 0
            call tytle
            newhead = .true.
         endif
c     Read cluster memberships from unit kc.
      if(normap) then
         read(kc)kk,lot,(r(i),i=1,nrank)
      else
         read(kc)kk,lot
      endif
      lit = lot -1
c     If this case is a member of the current cluster being printed,
c     store it in the row to be printed.
      if(lit .eq. lotpas) then
         ik = ik + 1
         is(ik) = k
         if(incase .gt. 0) then
c           Copy case ==> alph
            do 140 i=1,incase
               lchar = case(i,k)
  140       alph(i:i) = char
            id(ik) = alph
         endif
      endif
c     Is it time to print the row?
      if((ik .eq. m6) .or. ((k .eq. nsampl) .and. (ik .gt. 0) )) then
c        Print Cluster headers if new page or new cluster.
         if(newhead) then
            write(12,4050)lotpas
 4050       format(/,30x,'Members of Cluster', i3)
            write (12,3090)
            write(12,fmtx6)(idstrng,i=1,m6)
            write (12,3090)
            newhead = .false.
            lap = lap + 5
         endif
         if(incase .gt. 0) then
            write(12,fmtx7)(is(i),id(i),i=1,ik)
         else
            write(12,fmtx7)(is(i),i=1,ik)
         endif
         ik = 0
         lap = lap+1
      endif
  150 continue
  160 continue
      return
c
      end
c
      subroutine mapper(matrix,r,b)
      implicit real*8 (a-h,o-z)
      real*8 likly
      character*60 fmt
      logical normap
      common/a/likly,xlimit,nvbles,nsampl,ntypes,nsemi,nparam,normap,kx
      common/b/convrg,critrn,ratio,htest,oldlik,en,bigest,thresh,guard,
     1enmin,title(36),fmt,mxhier,maxtyp,maxav,nclust,kscrat,inp,kpun,kv,
     2 lot,lack,luck,iprint,ltime,maxt,iter,iterx,ltye,itmc,nrank,ivalue
     3,method,incase,ntyp(20)
c     The next three lines are necessary because PowerStation
c     FORTRAN cannot pass real arrays to character arrays in subr. call.
      logical*1 kblk,kper,kstr,lim,lookup(20),matrix(120,50),mlog
      character*1 mchar,blnk,clim,str
      equivalence (mlog,mchar),(kblk,blnk),(kstr,str),(lim,clim)
      dimension  r(nrank),b(21)
      data  lookup/'1','2','3','4','5','6','7','8','9','A','B','C','D',
     1 'E','F','G','H','I','J','K' /,kblk/' ' /,kper/'.' /,kstr/'*' /
c     This routine reads the file of discriminant scores and produces a
c     semi-graphical display on the printout.
      kc=3
      if(nrank.le.1) return
      ltypes=ntypes-1
      do 10 i=1,2
      i1=i+1
      if(i1.gt.nrank) go to 11
      do 10 j=i1,nrank
      do 2 m=1,48
      do 2 n=1,116
2     matrix(n,m)=kblk
      do 3 m=1,48
3     matrix(1,m)=kper
      do 4 n=1,116
4     matrix(n,47)=kper
      rewind kc
      do 5 k=1,nsampl
      read(kc)l,key,(r(ir),ir=1,nrank)
      m=5.0*r(i)+58.
      n=-3.0*r(j)+23.
      m= max0(2,min0(m,116))
      n= max0(2,min0(n, 46))
      lim=lookup(key-1)
      mlog = matrix(m,n)
c     if (matrix(m,n).ne. kblk) go to 5   [Not valid for PowerStation]
      if (mchar .eq. blnk) go to 5
c     if (matrix(m,n).ne. lim) lim=kstr   [Not valid for PowerStation]
      if (mchar .eq. clim) clim = str
5     matrix(m,n)=lim
      call tytle
      write (12,13) j,ltypes
      ly=8
      lb=1
      do 8 n=1,47
      lb=lb+1
      go to (7,7,6), lb
6     lb=0
      ly=ly-1
      write (12,14) ly,(matrix(m,n),m=1,116)
      go to 8
7     write (12,15) (matrix(m,n),m=1,116)
8     continue
      xd=-11.0
      do 9 m=1,21
      xd=xd+1.0
9     b(m)=xd
10    write (12,16) (b(m),m=1,21),i
11    continue
      return
c
12    format (6x,i4,10f7.2)
13    format (5x,'Discriminant Function',i2,50x,'Number of Types=',i2)
14    format (i4,116a1)
15    format (4x,116a1)
16    format (f12.0,20f5.0/45x,'Discriminant Function',i2)
      end


      function chisqp (chisq,ndfree)
      implicit real*8 (a-h,o-z)
      a=chisq*0.5
      t=dexp(-a)
      nn=(ndfree/2)*2
      if (nn-ndfree) 2,1,2
1     p=t
      t=a*t
      kw=4
      ak=2.0
      go to 3
2     c=dsqrt(chisq)
      p=2.0*area(c)
      t=0.7978846*c*t
      ak=1.5
      kw=3
3     if (ndfree-2) 6,6,4
4     do 5 i=kw,ndfree,2
      p=p+t
      t=(t*a)/ak
5     ak=1.0+ak
6     chisqp=p
      return
      end
      function area (w)
c     area above w on normal curve
      implicit real*8 (a-h,o-z)
      x=w*.70710678118654752440
      kot=2
      if (x) 1,1,2
1     kot=1
      x=-x
2     area=0.5/((((((.0000430638*x+.0002765672)*x+.0001520143)*x+.009270
     15272)*x+.0422820123)*x+.0705230784)*x+1.0)**16
      go to (3,4), kot
3     area=1.0-area
4     return
      end
      subroutine uplow(a,n,mode)
c     changes a symmetric matrix from upper to lower triangular form
c     or vice versa.
c     a = an n*n symmetric matrix stored in triangular form as a
c     one-dimensional array of n*(n+1)/2 elements.
c
c     mode = +1 if the input has ibm storage mode 1
c     mode = -1 if the input has the opposite half of ibm storage mode 1
c                                 examples
c              mode = +1                               mode = -1
c     a(1)   a(2)   a(4)   a( 7)              a(1)   a(2)   a(3)   a( 4)
c            a(3)   a(5)   a( 8)                     a(5)   a(6)   a( 7)
c                   a(6)   a( 9)                            a(8)   a( 9)
c                          a(10)                                   a(10)
c
      real*8 a,t
      dimension a(1)
      if(n.lt.3) go to 5
      node= (mode+3)/2
      nn=n*(n+1)/2
      n2=n/2
      nn2=nn/2
      go to (1,3),node
    1 k=nn+1
c     reverse the order of the variables by
c     interchanging a(i,j) with a(n+1-j,n+1-i)
      ij=0
      m=n
      do 6 i=1,n2
      m=m-1
      k=k-i
      l=k
      do 2 j=i,m
      ij=ij+1
      t=a(ij)
      a(ij)=a(l)
      a(l)=t
    2 l=l-j
    6 ij=ij+i
      go to (3,5),node
    3 l=nn
c     reverse a, interchanging a(i) with a(nn+1-i)
      do 4 i=1,nn2
      t=a(i)
      a(i)=a(l)
      a(l)=t
    4 l=l-1
      go to (5,1),node
    5 return
      end
      function squad (s,b,mx,ls,lb)
      implicit real*8 (a-h,o-z)
      dimension s(1), b(1)
c     this program calculates the quadratic form squad=bsb/2
c     where s is a semi-matrix in upper triangular form,b is a vector,
c     mx= order of matrices
c     the matrices start at s(ls+1) and b(lb+1)
      squad=0.0
      ij=ls
      n=mx+lb
      if (mx-1) 5,4,1
1     nf=lb+1
      na=n-1
      do 3 i=nf,na
      ij=ij+1
      cute=b(i)*s(ij)/2.0
      do 2 j=i,na
      ij=ij+1
2     cute=cute+s(ij)*b(j+1)
3     squad=squad+cute*b(i)
4     squad=squad+s(ij+1)*b(n)**2/2.0
5     return
      end
      subroutine syminv (symtrx,norder,determ,rdumy,sdumy,untryd,lbefor)
      implicit real*8 (a-h,o-z)
      logical untryd
c     acm algorithm 150 by h. rutishauser
c     this program replaces a symmetric matrix by its inverse.
c     symtrx is the symmetric matrix, stored in upper triangular form as
c     a one-dimensional array starting at symtrx(lbefor+1)  .
c     rdumy and sdumy are temporary storage vectors used in the program.
c     norder is the order of the matrix and the temporary storage areas.
c     determ= natural logarithm of the determinant
      dimension symtrx(1), rdumy(1), sdumy(1), untryd(1)
      determ=0.0
      do 1 i=1,norder
c     untryd(i)=1 before i has been used as a pivot, and =0 afterwards.
1     untryd(i)=.true.
      do 7 l=1,norder
c     search for pivot
      bigajj=0.0
      n=lbefor+1
      do 2 j=1,norder
      t=dabs(symtrx(n))
      if (.not.untryd(j).or.t.le.bigajj) go to 2
      bigajj=t
      k=j
      nk=n
2     n=n+norder-j+1
      if (bigajj.le.1.0e-50) symtrx(nk)=1.0e-50
      determ=determ+dlog(dabs(symtrx(nk)))
c      sweeps out the symmetric matrix symtrx for a selected pivot k,
c     resulting in the so-called partial inverse useful in stepwise
c     regression.
      untryd(k)=.false.
c     preparation of elimination step
      rdumy(k)=1.0/symtrx(nk)
      sdumy(k)=1.0
      symtrx(nk)=0.0
      n=lbefor+k
      if (k.eq.1) go to 4
      ka=k-1
      do 3 j=1,ka
      sdumy(j)=symtrx(n)
      rdumy(j)=symtrx(n)*rdumy(k)
      if (untryd(j)) rdumy(j)=-rdumy(j)
      symtrx(n)=0.0
3     n=n+norder-j
4     if (k.eq.norder) go to 6
      ka=k+1
      do 5 j=ka,norder
      n=n+1
      sdumy(j)=symtrx(n)
      rdumy(j)=-symtrx(n)*rdumy(k)
      if (.not.untryd(j)) sdumy(j)=-sdumy(j)
5     symtrx(n)=0.0
c     elimination proper
6     n=lbefor
      do 7 i=1,norder
      do 7 j=i,norder
      n=n+1
7     symtrx(n)=symtrx(n)+sdumy(i)*rdumy(j)
      return
      end
c
c     ..................................................................
c
c        subroutine eigen
c
c        purpose
c           compute eigenvalues and eigenvectors of a real symmetric
c           matrix
c
c        usage
c           call eigen(a,r,n,mv)
c
c        description of parameters
c           a - original matrix (symmetric), destroyed in computation.
c               resultant eigenvalues are developed in diagonal of
c               matrix a in descending order.
c           r - resultant matrix of eigenvectors (stored columnwise,
c               in same sequence as eigenvalues)
c           n - order of matrices a and r
c           mv- input code
c                   0   compute eigenvalues and eigenvectors
c                   1   compute eigenvalues only (r need not be
c                       dimensioned but must still appear in calling
c                       sequence)
c
c        remarks
c           original matrix a must be real symmetric (storage mode=1)
c           matrix a cannot be in the same location as matrix r
c
c        subroutines and function subprograms required
c           none
c
c        method
c           diagonalization method originated by jacobi and adapted
c           by von neumann for large computers as found in 'mathematical
c           methods for digital computers', edited by a. ralston and
c           h.s. wilf, john wiley and sons, new york, 1962, chapter 7
c
c     ..................................................................
c
      subroutine eigen(a,r,n,mv)
      dimension a(1),r(1)
c
c        ...............................................................
c
c        if a double precision version of this routine is desired, the
c        c in column 1 should be removed from the double precision
c        statement which follows.
c
      double precision a,r,anorm,anrmx,thr,x,y,sinx,sinx2,cosx,
     1                 cosx2,sincs,range
c
c        the c must also be removed from double precision statements
c        appearing in other routines used in conjunction with this
c        routine.
c
c        the double precision version of this subroutine must also
c        contain double precision fortran functions.  sqrt in statements
c        40, 68, 75, and 78 must be changed to dsqrt.  abs in statement
c        62 must be changed to dabs. the constant in statement 5 should
c        be changed to 1.0d-12.
c
c        ...............................................................
c
c        generate identity matrix
c
    5 range=1.0d-12
      if(mv-1) 10,25,10
   10 iq=-n
      do 20 j=1,n
      iq=iq+n
      do 20 i=1,n
      ij=iq+i
      r(ij)=0.0
      if(i-j) 20,15,20
   15 r(ij)=1.0
   20 continue
c
c        compute initial and final norms (anorm and anormx)
c
   25 anorm=0.0
      do 35 i=1,n
      do 35 j=i,n
      if(i-j) 30,35,30
   30 ia=i+(j*j-j)/2
      anorm=anorm+a(ia)*a(ia)
   35 continue
      if(anorm) 165,165,40
   40 anorm=1.414d0*dsqrt(anorm)
      anrmx=anorm*range/dble(n)
c
c        initialize indicators and compute threshold, thr
c
      ind=0
      thr=anorm
   45 thr=thr/dble(n)
   50 l=1
   55 m=l+1
c
c        compute sin and cos
c
   60 mq=(m*m-m)/2
      lq=(l*l-l)/2
      lm=l+mq
   62 if(dabs(a(lm))-thr) 130,65,65
   65 ind=1
      ll=l+lq
      mm=m+mq
      x=0.5*(a(ll)-a(mm))
   68 y=-a(lm)/dsqrt(a(lm)*a(lm)+x*x)
      if(x) 70,75,75
   70 y=-y
   75 sinx=y/dsqrt(2.0*(1.0+(dsqrt(1.0-y*y))))
      sinx2=sinx*sinx
   78 cosx=dsqrt(1.0-sinx2)
      cosx2=cosx*cosx
      sincs =sinx*cosx
c
c        rotate l and m columns
c
      ilq=n*(l-1)
      imq=n*(m-1)
      do 125 i=1,n
      iq=(i*i-i)/2
      if(i-l) 80,115,80
   80 if(i-m) 85,115,90
   85 im=i+mq
      go to 95
   90 im=m+iq
   95 if(i-l) 100,105,105
  100 il=i+lq
      go to 110
  105 il=l+iq
  110 x=a(il)*cosx-a(im)*sinx
      a(im)=a(il)*sinx+a(im)*cosx
      a(il)=x
  115 if(mv-1) 120,125,120
  120 ilr=ilq+i
      imr=imq+i
      x=r(ilr)*cosx-r(imr)*sinx
      r(imr)=r(ilr)*sinx+r(imr)*cosx
      r(ilr)=x
  125 continue
      x=2.0*a(lm)*sincs
      y=a(ll)*cosx2+a(mm)*sinx2-x
      x=a(ll)*sinx2+a(mm)*cosx2+x
      a(lm)=(a(ll)-a(mm))*sincs+a(lm)*(cosx2-sinx2)
      a(ll)=y
      a(mm)=x
c
c        tests for completion
c
c        test for m = last column
c
  130 if(m-n) 135,140,135
  135 m=m+1
      go to 60
c
c        test for l = second from last column
c
  140 if(l-(n-1)) 145,150,145
  145 l=l+1
      go to 55
  150 if(ind-1) 160,155,160
  155 ind=0
      go to 50
c
c        compare threshold with final norm
c
  160 if(thr-anrmx) 165,165,45
c
c        sort eigenvalues and eigenvectors
c
  165 iq=-n
      do 185 i=1,n
      iq=iq+n
      ll=i+(i*i-i)/2
      jq=n*(i-2)
      do 185 j=i,n
      jq=jq+n
      mm=j+(j*j-j)/2
      if(a(ll)-a(mm)) 170,185,185
  170 x=a(ll)
      a(ll)=a(mm)
      a(mm)=x
      if(mv-1) 175,185,175
  175 do 180 k=1,n
      ilr=iq+k
      imr=jq+k
      x=r(ilr)
      r(ilr)=r(imr)
  180 r(imr)=x
  185 continue
      return
      end
c
c     ..................................................................
c
c        subroutine nroot
c
c        purpose
c           compute eigenvalues and eigenvectors of a real nonsymmetric
c           matrix of the form b-inverse times a.  this subroutine is
c           normally called by subroutine canor in performing a
c           canonical correlation analysis.
c
c        usage
c           call nroot (m,a,b,xl,x)
c
c        description of parameters
c           m  - order of square matrices a, b, and x.
c           a  - input matrix (m x m).
c           b  - input matrix (m x m).
c           xl - output vector of length m containing eigenvalues of
c                b-inverse times a.
c           x  - output matrix (m x m) containing eigenvectors column-
c                wise.
c
c        remarks
c           none
c
c        subroutines and function subprograms required
c           eigen
c
c        method
c           refer to w. w. cooley and p. r. lohnes, 'multivariate pro-
c           cedures for the behavioral sciences', john wiley and sons,
c           1962, chapter 3.
c
c     ..................................................................
c
      subroutine nroot (m,a,b,xl,x)
      dimension a(1),b(1),xl(1),x(1)
c
c        ...............................................................
c
c        if a double precision version of this routine is desired, the
c        c in column 1 should be removed from the double precision
c        statement which follows.
c
      double precision a,b,xl,x,sumv
c
c        the c must also be removed from double precision statements
c        appearing in other routines used in conjunction with this
c        routine.
c
c        the double precision version of this subroutine must also
c        contain double precision fortran functions.  sqrt in statements
c        110 and 175 must be changed to dsqrt.  abs in statement 110
c        must be changed to dabs.
c
c        ...............................................................
c
c     compute eigenvalues and eigenvectors of b
c
      k=1
      do 100 j=2,m
      l=m*(j-1)
      do 100 i=1,j
      l=l+1
      k=k+1
  100 b(k)=b(l)
c
c        the matrix b is a real symmetric matrix.
c
      mv=0
      call eigen (b,x,m,mv)
c
c     form reciprocals of square root of eigenvalues.  the results
c     are premultiplied by the associated eigenvectors.
c
      l=0
      do 110 j=1,m
      l=l+j
  110 xl(j)=1.0/dsqrt(dabs(b(l)))
      k=0
      do 115 j=1,m
      do 115 i=1,m
      k=k+1
  115 b(k)=x(k)*xl(j)
c
c     form (b**(-1/2))prime * a * (b**(-1/2))
c
      do 120 i=1,m
      n2=0
      do 120 j=1,m
      n1=m*(i-1)
      l=m*(j-1)+i
      x(l)=0.0
      do 120 k=1,m
      n1=n1+1
      n2=n2+1
  120 x(l)=x(l)+b(n1)*a(n2)
      l=0
      do 130 j=1,m
      do 130 i=1,j
      n1=i-m
      n2=m*(j-1)
      l=l+1
      a(l)=0.0
      do 130 k=1,m
      n1=n1+m
      n2=n2+1
  130 a(l)=a(l)+x(n1)*b(n2)
c
c     compute eigenvalues and eigenvectors of a
c
      call eigen (a,x,m,mv)
      l=0
      do 140 i=1,m
      l=l+i
  140 xl(i)=a(l)
c
c     compute the normalized eigenvectors
c
      do 150 i=1,m
      n2=0
      do 150 j=1,m
      n1=i-m
      l=m*(j-1)+i
      a(l)=0.0
      do 150 k=1,m
      n1=n1+m
      n2=n2+1
  150 a(l)=a(l)+b(n1)*x(n2)
c
c     At this point, the routine differs from the standard ibm nroot.
c     The eigenvectors are not normalized to unit length. As a result,
c     instead of (x-prime)*(x)=i, we have (x-prime)*b*(x)=i , so that
c     the discriminant scores will have unit variances within-groups.
      n2 = m*m
      do 160 i=1,n2
  160 x(i)=a(i)
      return
      end